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Abstract 

We develop Polnarev's analytic method to calculate of the polarization power spectrum of the 
cosmic microwave background radiation (CMB) generated by cosmic relic gravitational waves. In 
this analytic approach the physics involved in this process is more transparent. Consequently, the 
effects due to various elements of physics can be isolated easily. In numerically calculating the evo- 
lution of the gravitational waves during the radiation-matter transition, the WKB approximation 
for the scale factor has been taken. To describe more precisely the decoupling process, we have in- 
troduced an analytic expression for the visibility function, consisting of two pieces of half-gaussian 
curves. We also include the damping effects on the power spectrum at small scale up to the second 
order of the tight coupling due to the collisions. An analytic polarization spectra C xx has been 
obtained with following several improvements over the previous results. 1. The approximate ana- 
lytic result is quite close to the numerical one evaluated from the cmbfast code, especially, for the 
first three peaks of the spectrum that are observable. By using the analytic exact solution of relic 
gravitational waves from the sudden-change approximation, we have demonstrated the dependence 
of C xx on the dark energy and the baryon. 2. Our analytic half-gaussian approximation of the 
visibility function fits analytically better than the usual Gaussian model, and its time integration 
yields a parameter-dependent damping factor. This improves the spectrum AC XX /C xx ~ 30% 
around the second and third peaks. 3. The second order of tight coupling limit reduces the am- 
plitude of Cf x by 58%, comparing with the first order. 4. The influence of the power spectrum 
index of relic GW is such that a larger value of n-j- produces higher polarization spectra C xx . 
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1. Introduction. 



Studies on the anisotropy and polarization of CMB have made great progress, yielding im- 
portant information of cosmology. Recently, Wilkinson Microwave Anisotropy Probe (WMAP) 
observational results on the CMB anisotropy and polarization power spectrums [1, 2, 3, 4, 5, 6] 
agree well with the prediction of inflation of a spatially flat Universe with the nearly scale- 
invariant and Gaussian spectrum of primordial adiabatic perturbations. Inflationary expansion 
can generate two type of perturbations of spacetime metric: one is the scalar type (den- 
sity) of perturbations [7, 8, 9], and the other is the tensorial type (relic gravitational waves) 
[10, 11, 12, 13, 14, 7, 9]. These two kinds of perturbations will enter the Boltzmann equation 
for photons and influence CMB during the decoupling. Their impact on the anisotropy and 
polarization of CMB are different, especially their respective contributions are not completely 
determined theoretically. The major contribution to CMB anisotropies and polarization are 
from the density perturbations [15, 16]. However, the tensorial contribution is also important, 
especially in long wave-length range. Moreover, the magnetic type of polarization of CMB can 
only be generated by the tensorial perturbations, and it thus provides another channel to detect 
the relic gravitational waves besides the direct detection of laser interferometer [17, 18, 19, 20]. 

The power spectra of CMB polarization can be calculated by numerical method [21, 22], 
which gives rather precise predictions. But the semi-analytic method is also very helpful in an- 
alyzing the underlying physics and in revealing the dependence on the cosmological parameters 
[23]. A common treatment uses the spherical harmonic functions to expand the Boltzmann 
equation into a hierarchical set of equations for the multipole moments, then solves each of 
them separately [21, 24, 25, 26]. The other treatment was first suggested by Polnarev [14], 
using a basis of polarization vectors to decompose the Boltzmann equation, ending up with 
only two equations for the two unknown distribution functions, ( and (3, standing for the 
anisotropy and the polarization, respectively. This treatment is simpler and has been further 
used [27, 28, 29, 30, 31, 32, 33]. In this paper, we study the CMB polarization caused by 
the relic gravitational waves in the Polnarev framework. Our result is an analytic formula 
of the polarization power spectrum, which depends explicitly on the visibility function and 
on the spectrum of the relic gravitational waves at the decoupling. In our treatment for the 
ionization history through recombination we use a half-gaussian visibility function, which is 
more precise than the usual gaussian function. For the relic gravitational waves we use a WKB 
approximation. When doing integration for the Boltzmann equations, we find there are two 
kind of damping mechanism on the power spectrum: one is caused by the visibility function, 
which causes damping at small scale, and another is the second order tight coupling limit 
when solving the Boltzmann equations, which is only a simple factor damping in all scale. The 
latter is different from the case of the scalar perturbation sources, where the so-called "silk 
damping" , which only causes the damping for the power spectrum at the small scale. We also 
find that the finial power spectra depend sensitively on the state of tensor perturbations at the 
decoupling time (the ratio of the positive and negative modes) , which is complex and prevents 
us from getting an exact formula for the power spectrums. 

The organization of this paper is as follows: in Section 2 and 3, we write down the evolution 
equations of the polarization distribution function and the exact relations of Cf G and Cf c 
with it. For the following discussion, we calculate the evolution of the gravitational waves and 
build the visibility function models in Section 4 and 5. In Section 6, we have an approximately 
analytic calculation of the Boltzmann equations, and give the expression of fy, which follows 
the expression of the power spectrum functions Cf G and Cf°. At last, we give a conclusion 
and discussion in Section 7, where we main discuss the effect of tensor-scalar ratio r, the baryon 
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density ,dark energy density and the reionization process on the power spectrums. 



2. CMB Polarization. 

The polarized distribution function of photons is generally represented by a column vector 
/ = (It, I r ,U,V), and its components are related to the Stokes parameters: I — Ii + I r and 
Q = l\ — I r . An important property of the Stokes parameters is that, under a rotation 5 about 
the axis of propagation, the total intensity I and the parameter V are invariant, but Q and U 
transform as [34] 

( Q' \ _( cos25 sin25 \(Q\ 
\ U' ) ~ \ -sin25 cos 25 ) \ U J ' 

So (Q, U) together form a spin-2 field according to the coordinate transformation, and can be 
conveniently described by a 2 x 2 polarization tensor. For actual detections the photon come 
from the full sky of 2-sphere, which in the spherical coordinates (9, 0) has a metric 




9ab~ [ n ^;„2 n ) • ( 1 ) 

the polarization tensor is [30] 




Q(n) -U(n)sm9 
-U(n)sm9 Q{n)s\n 2 9 



2a I, (2) 



satisfying P ah = P ha , and g ab P ab = 0. 

During the era prior to the decoupling in the early Universe, the Thompson scattering of 
anisotropic radiation by free electrons can give rise to the linear polarization only, and does not 
generate the circular polarization V, so we only consider / = (//, I r , U). For the trivial case of a 
homogeneous and isotropic unpolarized radiation, the distribution is simply / = f (u)(l, 1, 0), 
where fo{v) = -t^jvfzi is the usual blackbody distribution function with temperature T. The 
combined effects of the Thompson scattering and the metric perturbations will yield linear 
polarizations of photons. The time evolution of the photon distribution function is determined 
by the equation of radiative transfer, essentially the Boltzmann equation [34], 

?L + h i — = -—®L - a -j) 0) 

dr) dx l drj dv 

where n l is the unit vector in the direction (9, 0) of photon propagation, q = UTn e a is the 
differential optical depth and has the meaning of scattering rate, ot = 6.65 x lCT 25 cm~ 2 is the 
Thompson cross-section, n e is the number density of the free electron, and 

1 rl rl-K 

3 = 47 / i dfI ' 'I W'PM^^tVkn,**,"^,*'), (4) 

where \i = cos 9, ji' = cos 9' and 

/ / uV 2 cos2(0' -0) -/i 2 cos 2(0' -(f)) //V sin 2(0' -0^ 
P= -/i 2 cos 2(0' -0) cos 2(0' -0) -/isin2(0'-0) j (5) 
v -2/V 2 sin 2(0' - 0) 2fi sin 2(0' - 0) 2/i// cos 2(0' - 0) 
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is the phase-matrix. The scattering term q(f—J) in Eq.(3) describes the effect of the Thompson 
scattering by free electrons, and the term — reflects the effect of variation of frequency 
due to the metric perturbations through the Sachs- Wolfe formula [35] 



\dv 
v drj 



ldhi 



2 drj 



—n l h j . 



In the presence of perturbations hij, either scalar or tensorial, the distribution function will 
be perturbed and can be generally written as 



/(M) = /o 




(6) 



where fi represents the perturbed portion. 

The perturbed flat Robertson- Walker (FRW) metric is 

ds 2 = a 2 (r]) drj 2 - (Sij + h ij )dx i dx j . (7) 

where r\ = J(ao/a)dt is the conformal time, and are the perturbations with \hij\ <C 1. 
In our context, we consider only the tensorial type perturbations h^, representing the relic 
gravitational waves. So they are symmetric = hji, traceless hu = 0, and transverse hijj = 0. 
Therefore, there are only two independent modes, corresponding to the + and x gravitational- 
wave polarizations. 



h + 4 + h>< 4- 



Taking the direction of propagation of the GW in the direction of z, k — z, then the polarization 
tensors for the GW satisfy 



e^rnrij 



sin ft cos. 



efjhihj = sin 2 9 sin 20. 



In cosmological context, it is usually assumed that the two components h + and h x have the 
same magnitude and are of the same statistical properties. To simplify the Boltzmann equation 
(3), for the = h + efj polarization, one writes the perturbed distribution function fi in the 
form [14] 

> / 1 \ a ( (l+/i 2 )cos20 

[oj 2 



For the h. 
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2 (l-y 2 ) cos 20 
h x e*j polarization, one writes fi in the form 



V 



1 + /i 2 ) cos 20 
4/i sin 20 



(8) 



/ 1 = |(l-/i 2 ) sin 20 



( ! ) 


< 


V o y 


\ 



-(l + /i 2 )sin 20 
— 4/x cos 20 



(9) 



where ( represents the anisotropy of photon distribution since ( oc Ii +I r = I, and f3 represents 
the polarization of photons since (5 oc l\ — l r = Q. Both ( and (5 are to be determined by solving 
the Boltzmann equation. For the = h + e^ polarization, one substitutes / into Eq.(3). Upon 
taking Fourier transformation, retaining only the terms linear in the perturbation h^, and 
performing the integration over d/i, one arrives at a set of two differential equations [14, 32, 33], 

rfln/o 



£ fc + [ikfi + q] Cfc 







(10) 
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Pk + [ikn + q](3 k = ^qj dfi' [(l + /i' 2 )* (3 k t 



;n) 



For the hij = h x e^ polarization, the resulting equations are the same as Eqs.(lO) and (11) 
with h + being replaced by h x . In Eq.(lO) and (11) = ( k + Pk, k is the wavenumber, 
Pk, and h\ are the Fourier modes of £, /3, and h + , respectively, and the over dot " • " denotes 
d/drj. In the following, for simplicity, we will omit the sub-index k of the function (5 k and 
/i fc . Moreover, we also drop the GW polarization notation, + or x, since both h + and h x are 
similar in computations. From the structure of Eqs.(lO) and (11), one can see that the h of 
GW in Eq.(lO) plays the role of a source for the anisotropies £, which in turn plays the role of 
a source for the polarization (3 in Eq.(ll). Our work is to find the solution of P, then calculate 
the CMB polarization power spectra. 

In the simple case of the long-wavelength limit with k — 0, this set of equations reduces to 

{ + s£ = h, (12) 
* + > = -:*>*■ (13) 

The solutions P and £ will be independent of the angle /i. In general case of k ^ 0, the function 
P and £ to be determined by Eqs.(10) and (11) will depend on //. As the right hand side of 
Eq.(ll) contains an integral over d/j,', it is difficult to give an exact solution. And one may 
expand ft and £ in terms of the Legendre functions 

^)=D 2I + 1 )WW- 
i 

i 

with the Legendre components 

Mv) = l£ i d t Mt(T,, ti )Pi(jM), (14) 

Pi(v) = lj' i dfiP(v,^)Pi^)- (15) 

The differential equations Eqs.(10) and (11) for £(r},[i) and /3(ry,/x) then become an infinite set 
of coupled differential equations for £i(r)) and Piirf). 



3. Decomposition of Polarization f3 into Electric and Magnetic Types 

From Eq.(8) for the definition of ( and P for the h + GW polarization , it is seen that the 
Stokes parameters can be be written as the following [30, 32], 

Q(M) = ^E( 2/ + 1 ) P Kcos^)(l + cos 2 ^)cos20 A; (16) 

C/(0,0) = ^5Z(2Z + l)P,(cos0)2cos0sin20 A- (17) 

One commonly uses the spherical harmonic functions V(j m ) as the complete orthonormal basis 
for scalar functions defined on the 2-sphere, such as the temperature anisotropies AT. For the 
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2x2 tensors defined on the 2-sphere, such as P ab in Eq.(2), the following complete orthonormal 
set of tensor spherical harmonics can be employed [30]: 



Y, 



G 

(lm)ab 



c I • 



(lm)ab 2 



Y(j m y_ ac € b -\- Y(i m y bc € a^j j 

where ":" denotes covariant derivative on the 2-sphere, N t = y / 2(Z-2)!/(/ + 2)!, and 



sin 6 
— l/siii0 



They satisfy 



/ 



dnY^ ab (h) Yy,^(h) = 0. 



(18) 
(19) 

(20) 

(21) 
(22) 



By construction one sees that Y^ m ^ ab is the gradients (electric type) and Y^ m ^ ab is the curls 
(magnetic type) of the ordinary spherical harmonics, and they represent the electric type and 
magnetic type components of the polarization, respectively. The polarization tensor can be 
expanded in terms of this basis as: 



oo I 

-o z=2m= _/ 

The expansion coefficients are given by 



Pg b (n) 



a (Zm)^n)ab' 



™) + a (lm)Y(lm)ab(n) 



a (lm) 



l -f dnP ab {n)Y^*{n) 



l (lm) 



dhP ab {n)Y { 



(23) 



(24) 



and calculation yields [30, 32]: 1 

a?m= \(Sm2 + S m ,-2)y/27T(2l + l) 

'(Z + 2)(Z + l)A-2 , 6(Z-1)(Z + 2)A 



X 



l(l - l)A+2 



(2/-l)(2/ + l) (2/ + 3)(2/-l) (2/ + 3)(2/ + l) 



2tt 



-(S m2 - 6 m> - 2 )[(l + 2)A-i + (/ - l)A+i]- 



(25) 



(26) 



lm 4 V (2Z + 1) ' 

where the polarization /3 shows up explicitly in the G and C components. Then the magnetic 
type of power spectrum Cp G is 2 

1 



GG 



2/4 

7T 

16 



(Z + 2)(Z + l)A- 2 + 6(Z-l)(Z + 2)A 



z(z - i)A+2 



(2Z-1)(2/ + 1) (2Z + 3)(2/-l) (2/ + 3)(2/ + l) 



(27) 



x There is a small mistake in the formulae (4.39), (4.40), (4.41) and (4.42) in Ref[30] and formulae (111), (114) and 
(115) in Ref[32], the coefficient of (2i+3)(2i-i) should be replaced by (2i+3)(2i~-i) ■ 

2 In the Ref.[30, 32], the bracket of [ ] should be replaced by the absolute value sign | |, which is for the values of /3 
may being imaginary numbers. 
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and similarly for C, cc . Note that these expressions are for the gravitational wave in the z 
direction with '+' polarization. Summing over all Fourier modes, and over both polarization 
states, one has the final result 



^GG = _J_ f 
1 167T J 



(Z + 2)(Z + l)A- 2 , 6(Z-1)(Z + 2)A , Z(Z-l)A +2 



(21 - 1)(2Z + 1) 
r cc _ 1 f 



(2Z + 3)(2Z-1) (2Z + 3)(2Z + 1) 

2 



(z + 2)A-i + (z-i)A+i 



2/ + 1 



2Z + 1 



fc 2 dife. 



Note that the cross-correlation power spectrum vanishes 

m=l G* C 
/^GC _ \ -> a lm a lm _ n 

6/ " 4-, 27TT" 0, 



(28) 
(29) 

(30) 



smce a (L) K (<W + 5m -2), while ag m) OC (6 m ,2 - S m - 2 ). 



4. Evolution of Gravitational Waves 

Now we consider the evolution of the relic GW, which is the source term of the CMB 
polarization in Eq.(10). For both polarizations, +, x, the equation of motion for the relic GW 
of mode k is the following: 

h + 2-h + k 2 h = 0, (31) 
a 

and the initial condition is taken to be 

h(n = 0) = h(k), h(r] = 0) = 0, (32) 

with 

— \h(k)f = P h (k) = A T (j-j , (33) 

where Ph{k) is the the primordial power spectrum of GW, At is the amplitude, k = 0.05 
Mpc -1 is the pivot wavenumber, and ut is the the tensor spectrum index. Inflationary models 
generically predicts tit ~ 0, a nearly scale-invariant spectrum. Later we will also see the influ- 
ence of n T on the CMB polarization. We have ignored a suppressing effect on the gravitational 
waves by the neutrinos free streaming [36, 37], which can slightly bring down the hight of the 
peak at small scales [26]. 

The equation (31) depends on the scale factor a(r)), which is determined by the Friedmann 
equation 

a 2 = Hi [tt r + aQ m + a 4 fi A ] , (34) 

where H is the present Hubble parameter, Q r , Q m , and are the present fractional densities 
for the radiation (including the photon and neutrino), matter, and dark energy, respectively. 
Given a set of these densities, one solves Eq.(34) for a(rf) numerically. Taking Q r = 8.36 x 
10~ 5 , Vt m = tt b + Q dm = 0.044 + 0.226, tt A = 0.73, we solve Eq.(34), and substitute the 
resulting a(rj) and a{rj) into Eq.(31), then the numerical solutions h(rj) and h[rj) are obtained 
straightforwardly. The resulting h(r)d) and h(r)d) at the decoupling time rjd are given as function 
of k in Fig.l and Fig.2, respectively. 
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Besides the numerical solution, we may use the simple approximations of a(n) transiting 
between different stages suddenly, 

!a r n, T) < r) e (radiation dominant), 

a m n 2 , n e < n < r\ E [matter dominant), (35) 
ain~ l , n > t]e (A dominant) 

where a r , a m , and a\ are constant numbers, and can be determined by jointing a(n) at n e and 
n E . In the ACDM model 3 with Q b = 0.044, Q dm = 0.226, Q A = 0.73, taking the redshift 
z e = 3234 yields the conformal time n e /riQ = 0.007 (where r/ is the present conformal time), 
and taking z E = 0.39 yields n E /rj = 0.894, and 7^/770 — 0.0195 is the decoupling time at the 
redshift z d = 1089. Eq.(31) has the analytic solution [38, 25, 39]: 

h(n) = Aojo(kn), (n < rj e ), (36) 

h{n) = A^nJn^A^kn) + A 2yi {kn)}, {n e < n < n E ), (37) 

with the coefficient 

A = (A T k nT ~ 3 ) 1/2 , (38) 
determined by the primordial power spectrum, and 

3kn e - kn e cos(2kn e ) + 2 sm(2kr) e ) 
Al ~ 2k^l ' (39) 

2- 2k 2 n 2 e - 2 cos{2kn e ) - kn e sin(2%) 
M ~ 2k^l • (40) 

These h(nd) and h(r]d) are plotted with solid lines in Fig.l and Fig.2, where the initial normal- 
ization h(k) = 1 has been taken. The figures show that the simple approximation is good only 
in long wavelengths, but it differs from the numerical one considerably in short wavelength 
(k > 10). 

As an improvement to (35), one can use the WKB approximation for a(n) [40, 26]. Since 
we are only interested in the gravitational waves at time r\d, and the dark energy is small and 
can be omitted at that time, so the scale factor can be approximated by 

a(r)=a e r(r + 2), (41) 

where r = (y/2 — l)n/n e , and a e is determined by a /a e = 1 + z e . When r <C 2, a(r) — > r, 
the radiation dominated stage, and when r ^> 2, a — > r 2 , the matter dominated stage. a(r) 
transits between these two stages is smooth. Then the evolution of the gravitational waves 
become 

h" + 2-ti + r 2 h = 0, (42) 
a 

where r = kn e /(\/2 — 1), the prime denotes d/dr. This equation has an analytic solution 
discussed in Ref. [26], which is rather complex. In this paper, we employ Eq.(41) as an approx- 
imation of a(n) to find the numerical solution of Eq.(42). The resulting h(n d ) and h(n d ) in this 
WKB approximation are plotted with the dashed lines in Fig.l and Fig.2, which show that 
the results are very good when comparing with the numerical ones, and the difference of them 
1%. The approximation of (41) is simpler than the numerical a(n), and has better precision 
than the simple approximation of (35), and will be used to calculate the CMB polarization 
power spectrum. 



3 In all this paper, we choose the present Hubble parameter ho = 0.72. 
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5. The Model of Decoupling History. 



Consider the decoupling history of the Universe. Before the decoupling, the ionized baryons 
are tightly coupled to photons by Thompson scattering. Once the temperature falls below a 
few eV, it becomes favorable for electrons and ions to recombine to form neutral atoms. As the 
number of charged particles falls, the mean free path of any given photon increases. Eventually, 
the mean free path becomes comparable to the horizon size and the photon and baryon fluids 
are essentially decoupled, and the CMB photons last scatter. One can solve the ionization 
equations during the recombination stage to obtain the visibility function V(r)), which describes 
the probability that a given photon last scattered from a particular time [41, 42], depending 
on the cosmological parameters, especially and the present Hubble parameter Hq [23]. In 
terms of the optical depth k, this visibility function is given by 

V{ V ) = q{r])e- K ^°' v \ (43) 

satisfying 

rvo 

v(n)d n = i, (44) 



f 

JO 
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where the optical depth function K,(r)o,r)) is related to the differential optical depth q(rj) by 
q{rj) = —dK(r)o,r))/dr). Fig.3 shows the profile of V(r)) from the numerical result by cmbfast, 
which is sharply peaked around the last scattering. In practical calculation it is usually fitted 
by a narrow Gaussian form [24, 26] 

V( V ) = V( m )e X p^-^^j, (45) 

where V(r)d) is the amplitude at the the decoupling time rjd, and Arj d is the thickness of 
decoupling. The analysis of the WMAP data [3] gives the redshift thickness of the decoupling 
Az d = 195 ± 2, which corresponds to Ar] d /r]Q = 0.00143. Then, taking V(r] d )r]o — 279 and in 
(45) yields a fitting shown in Fig.3, which has large error on both side of 77^, compared with 
the numerical one. To improve the fitting of V(i]), we take the following analytic expressions, 
which consists of two half-gaussian functions, 



( ?7 ~ Vd Y 



^)=n%)exp(-^^], (77<77 d ); (46) 



V(v) = V(v d )ew{-%^-y (v>Vd); (47) 

with Ar] dl /r] = 0.00110, Ar) d2 /r} = 0.00176, and (Arj d i + Ar] d2 )/2 = Ar/ d , satisfying the 
constraint of (44). Fig.3 shows that the half-gaussian model fits the numerical one much 
better than the Gaussian fitting, especially, on the left-hand side of the peak the area enclosed 
between the curves of the Gaussian and the half-gaussian is about ~ 11% of the total area 
enclosed under the curve V{rj). As shall be seen later, this difference in V(rj) will subsequently 
cause a variation in the hight of the polarization spectrum. The expressions (46) and (47) will 
be used to calculate the approximate analytic polarization power spectrum, which turns out 
to depend more sensitively on the smaller time interval Ar] d i- 



6. Analytic Solution. 
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We start to look for the analytic solution of Eqs.(lO) and (11). Since the blackbody spec- 
trum f(vo) in the Rayleigh- Jeans zone has the property dl ° 1 ^°- > 1, these equations reduce 
to 

£ + [ikfi + q]i = h, 
+ [ ikl x + q](3 = ^ £ dp! 



;i+// 2 ) 2 /3~(i-/m 



(49) 



In Eq.(48) the GW h play the role of source for the anisotropies £, while the term g£ causes £ 
to damp. The formal solution of Eq.(48) is 

^{r})= r hi^e-^'V^'-^dr]'. (50) 
j o 

In Eq.(49) the integration over /i contains the functions (5 and £. Using ^ and A defined in 
Eqs.(15) and (14), then Eq.(49) is 

$ + [ikfi + q]P = qG, (51) 

where 

G(v) = |a + \a, + 1 A - ^e. + i& - -56.. 

One might write down a formal solution 

/Jfo) = T G{r]')q(r]')e- K ^'V k ^'-^dr]', (52) 

and set the the time 77 in the above to be the present time 770, 

(3( Vo ) = / G{rf)V(rf)e ik ^^drf, (53) 
Jo 

where V(rj') = q< K r]')e~ K ^ ri0 ' ri '^ is the visibility function. However, the difficulty with the inte- 
gration (53) is that the integrand G contains A and ^ up to I — 4, which are not known 
yet. 

One uses the Legendre expansion and write Eqs.(48) and (49) as the following hierarchical 
set of equations: 

io = -q^o ~ iKi + h, (54) 

3 { 3 5 3 1 1 \ 

A) = - ik(3i + q (—/3 4 + -fa ~ ^£4 + ^2 - — £0) , (55) 

6 = -?6 - + (' + ^i+d » /° r z ^ !» ( 56 ) 

A = -<zA - ^7 [/A-i + (i + i)A+i] , /«• i> 1. (57) 

In the tight coupling limit with q — > oo, the equations reduce to 

£o + g£o = (58) 

A) + y^A) = -7^0, (59) 
6 = A = o, z > l. (60) 
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Then the source function G{rj) reduces to G — (7(3 — £ )/10, and satisfies the equation: 
and the formal solution is 

Glrj) = -1- [ V h{rf)e-& K W>drf. (62) 
10 Jo 

Substitute this expression of G into Eq.(53), yields the formal solution for the polarization in 
the tight coupling limit: 



Pirn) = r V M) f h( V ")e-^',v") dr] » 

Jo ID Jo 



e ikKv'-Vo) dr] > 



= _J_ f m dr]'V(ri')e ik M-™) H dr/'h(ri")e-& K W' )+ & K W\ (63) 
10 jo jo 

where k(t)) = K(r)o,r)). 

Note that this approximate result of the tight coupling applies only on scales much larger 
than the mean free path of photons. On smaller scales the effects of photon diffusion will cause 
some damping in the anisotropy and polarization. To take care of this effect, we now expand 
Eqs.(54)-(57) to the second order of the small parameter 1/q << 1, which has the meaning of 
the mean free path, and arrive at, 

io = -q£o - ikfr + h, (64) 

ii = -?6 - y£o, (65) 

6 = 0, I > 2. (66) 

Putting £ oc e < / w *' and & oc e*/^" and substituting into Eqs.(64) and (65), ignoring vari- 
ations of uj on the expansion scale a/a, neglecting h which is nearly zero at low frequency, 
shown in Fig. 2, one gets 

« = ±^ + » 

Thus £o will acquire an extra damping factor e~ I qdv y independent on the wavenumber k. This 
feature is different from the case of the scalar perturbations, where the damping is strong on 
the small scales [43]. For the polarization /3, we only keep the tight coupling limit with the 
equation, 

Po - -y^o ~ ^o, 

which follows that /3 also gets a damping of e~ -I qdri . Thus, taking into the account of the 
diffusion effect, G in (62) will acquire an extra damping factor exp(— k(t))), and (63) becomes 

p( Vo ) = --L f m drfV(rf)e ik M-*>) f' drf'h(rf')e-^ K ^-^ K ^\ (67) 
10 jo jo 

In the above the function exp (— ^k(t}")) ~ for rj < rj d , and exp (— ^n(rj")) ~ 1 for 
i] > r/d, so it can be approximated as a step function exp (— ^k(t)")) ~ &(r)" — rjd), and 
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moreover, the visibility function V(r)') is also peaked about the decoupling r\ d . Therefore, as 
an approximation, one can pull the h(r)") out of the integration / drj", 

P( ! n Q ) = -L r dr]V(ri)e lk,l ^ m) h(ri) f dr/V^ " 1 ^. (68) 
10 Jo Jo 

Define the integration variable x = K,(r]') / K,(rj) to replace the variable 77'. For V(rj) is peaked 
around the r] d with width Arj d , one can take drj' = ^-Ar] d as an approximation, then 

/3( m ) = -^-Arj d [ V0 drjVMe^-^hirj) [°° — e -^ x e~^\ (69) 
10 Jo Ji x 

Substituting this into Eq.(15)and using the expansion formula 



Akr cos 



= J E(2l + l)i l j l (kr)P l (cos6), 

1=0 



one arrives at the expression for the component of the polarization 

I /"^o r°° dx 3 7 

PiiVo) = -TzAVdi 1 / drjV^h^jtik^-rto)) _ e -io«(^ e -io«W. (70) 
10 Jo J\ x 

Now look at the integration / drj involving V(r)), which has a factor of the form e~ a( - v ~ Vd ^ . As 
a stochastic quantity, the time- derivative of GW h{rj) contains generally a mixture of oscillating 
modes, such as e %kri and e~ tkr >, and so does the spherical Bessel function ji(k(r) — 770)). Thus 
h(r))ji(k(r) — ?7o)) generally contains terms oc e~ ibk ^ Vo \ where b G [—2,2]. Using the formula 

r e~ ay2 e lbky dy = e~^r f°° e - ay2 dy, 

J —oo J —oo 

the integration is approximately 

"* d V V( V )h( v ) 3l (k( v - 770)) ~ e-^ kA ^ 2 h( Vd ) 3l (k( Vd - 770)) r d v V( V ). (71) 

Jo 

When the half-gaussian visibility function V(rj) of Eq.(46) and (47) is used, the integration is 
approximated by 

P dvVWhWMHr, - 770)) » ^ [ e ^( fcA ^) 2 + e -(^ d2 ) 2 j h(rj d )ji(k(r] d - rj )) P ^(77). 

JO Z L J JO 

(72) 

In the above a can take values in [0, 2], depending on the phase of h(r))ji(k(r) — 770)). Here we 
will take a as a parameter. This gives another damping depending on both Ai] d and k. This 
damping of the anisotropics £ and (3 of CMB is caused by Thompson scattering of photons by 
free electrons. During the recombination around the last scattering, the visibility function V 
is narrowly centered around the the time 77^ with a width Ar/ d , so the smoothing by Thompson 
scattering is effectively limited within the interval Arj d . Thus a wave of anisotropies e lkv will 
be damped in this interval by an factor e - ( fcAr/d ) . The longer the interval Ar/ d is, the more 
damping the wave suffers. In fact Ar] d can be viewed as the thickness of the last scattering 
surface. Those waves with wavelength A shorter than Ar/ d will be effectively damped by a 
factor e ~( 2wA ^/ A ) 2 , the same as before. Thus the shorter the wave length is, the more damping 
for the wave. 
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The remaining integrations in (3i is 

i m drjV(rj) r ^e-^*e^ = f f = -!*-■ (73) 

J o Ji x Jo Ji x 17 3 

We like to point out that this amplitude is the outcome from the second order of the tight- 
coupling limit, while the first order of the tight-coupling limit with G being given in (62) would 
yield a result f In f in Ref.[26]. Thus finally (70) yields 

1 20 

= -t= In -i l ^ d h{m)3mr, d - r)o))D(k), (74) 



D(k) = - [ e -«(^ d i) 2 + e -a(kA Vd2 f} _ ^ 
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where 

2 

for the half-gaussian visibility function. For the gaussian visibility function one would have 
D(k) = e - a( - kA ^ 2 . 

Substituting this back into Eqs.(28) and (29) yields the polarization spectra 
1 / 1 90\ 2 r 

Cf X = — (^ ln y) J P 2 xMv d -Vo))\h(v d )\ 2 ^V 2 d D\k)edk, (76) 

where "X" denotes "G" or "C", the type of of the CMB polarization, for the electric type 

n . , (/ + 2)(/ + l) . . . 6(/-l)(/ + 2) . . . 1(1-1) . .. 

PGI{X) = (21-I)(2l + If - 2{X) ~ (21-1)(21 + 3) MX) + (21 + 3)(21 + 1) JI+2{X) ' (7?) 

and for the magnetic type 

Pci{x) = - (78) 

The result (76) is similar to the result in Ref.[26] (realizing that our (2n) 3 Fourier conven- 
tions differ from theirs) if we identify Cf G = Cei/2 and Cf c = Cbi/2, but here the coefficient 
is Yf m y> smaller than that in Ref. [26] since we have also included the diffusion effect on the 
source G. So the the second order of the tight-coupling limit reduces the amplitude by about 
~ 58%. Another difference is the damping factor D{k). Besides, the parameter a is taken to 
be in the range [0, 2], instead of a fixed a = 1/2 as in Ref. [26]. 

To completely determine Cf x , we need to fix the normalization of the initial amplitude 
h(rid) in Eq.(76). What has been observed is the CMB temperature anisotropies, which gen- 
erally has contributions from both the scalar and tensor perturbations. The ratio of the 
contributions 

>-m 

has not been fixed observationally, and only constraints can be given. Based upon the obser- 
vations of Ly-a forest power spectrum from the SDSS, 3-year WMAP, supernovae and galaxy 
clustering, one can give a constraint of r < 0.22 at 95% C.L. (< 0.37, at 99.9% C.L.) [44]. We 
take it as a parameter. WMAP observation [4] gives the scalar perturbation power spectrum 

P R (k ) =2.95x 10- g A(k ), (80) 

with k = O.OSMpc^ 1 and the amplitude A(k ) = 0.8. Taking the scale-invariant spectrum 
with fix = in (33), then the amplitude A T in (33) depends on r. For instance, if r = 1 is 
taken, then A T = 2.46 x 10~ 9 , and smaller r will yield smaller A T accordingly. 
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7. Results and Discussions 



Damping Effects due to Visibility Function: 

The power spectra of Cf G and C[ ;c , calculated from our analytic formulae (76) and from the 
numerical cmbfast, have been shown in Figs. 4 and 5, respectively. The approximate analytic 
result is quite close to the numerical one evaluated from the cmbfast code, especially, for the 
first three peaks of the spectrum that are observable. One sees that Cf G and Cf c at large I 
sensitively depend on the visibility function V, that is, on the factor D{k). In particular for 
the electric polarization spectrum, our half-gaussian model D{k) in (75) with a = 1.7 gives 
very good fitting with the third peak being very close to the numerical one, and the second 
peak being higher. For a = 2 the second peaks of the spectra are good, but the third peaks 
are a bit too low. On the other hand, the Gaussian model D(k) = e~ a( - kAlld ^ with a = 2 yields 
a power spectrum too low. The reason is for Arj d > Arjdi- The larger value of a, the larger 
damping of the power spectrum on the small scale. 

The height of the power spectrums: 

Eq. (76) shows that the height of C xx depends on the amplitude \h(rj d )\ at the decoupling 
time. As has been discussed earlier, for the scale-invariant power spectrum nx — 0, \h(r) d )\ is 
directly related to the tensor-scalar ratio r in Eq.(79). A larger r yields a larger \h(n d )\ and 
a higher polarization. Currently the observations have only given an upper limit of r < 0.22 
(95% CL.) [44]. In Fig. 6 we have plotted the analytic formula (76) for Cf , for three values 
r = 0.3, 0.1, 0.01, respectively, whereby also plotted are the one-sigma sensitivity estimates of 
the near-term projects, WMAP and Planck satellites [45, 46, 47]. The WMAP estimates are 
based on measured noise properties of the instrument for an 8- year of operation, and the Planck 
estimates are based on noise measurements from the test-bed High Frequency Instrument for 
a 1.2-year of operation. Fig. 6 clearly shows that the magnetic polarization for the models with 
r > 0.1 can be detected by the Planck, but difficult for the WMAP. 

However, in this discussion, we have not considered the effect of cosmic reionization pro- 
cess, which is directly related to the galaxy formation. The recent WMAP result [5] tends to 
give the optical depth of reionization k t = 0.09 ± 0.03. Thus the visibility function V(n) will 
have another peak around the late time 77/770 ~ 0.27 besides the cosmic decoupling, and will 
give an extra contribution to Cf c , correspondingly. At present, the reionization process is not 
well understood yet, and it is difficult to give an analytic formula for this process. Using the 
numerical cmbfast including the reionization effect, we have plotted Cf in Fig. 7, where an 
extra peak of Cf c at I ~ 6 is seen. On the observational side, a number of other projects are 
currently underway, such as CBI[48], DASI[49], CAPMAP[50], BOOMERANG[51], emphasiz- 
ing on the CMB magnetic polarization. The future projects CLOVER[52] and QUIET[53] are 
expected to detect the magnetic polarization for r > 0.01, and the project CMBPOL [54]even 
for r > 10~ 3 . 

The influence of width of decoupling: 
Besides the gravitational waves, C xx also directly depend on the thickness of the decoupling 
An d , and on the damping function D(k). A smaller An d makes the power spectrum having 
smaller height. This effect is obvious at the large scale (small I). At the small scale (large 
I), the effect is complicate since both Ar/ d and the damping factor D(k) will influence the 
spectrum. The expression of D(k) shows that, for a fixed k, the smaller An d leads to larger 
D(k). The thickness of the decoupling is mainly determined by the baryon density Qt, of the 
Universe. In the flat ACDM Universe, increasing fl b will slightly enhance the decoupling speed, 
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which will make Arja becoming smaller [42]. For example, a fitting formula can be used for 
the optical depth function in the ACDM Universe [23] 

^o,z) = ^ 1 ( T ^) C2 , 800 <^< 1200, (81) 

where c\ = 0.43 and C2 = 16 + 1.8 In fife. This function is only dependent on fife. The visibility 
function V(rf) = ^e~ K is peaked around at r/ d . A larger fi b corresponds to a narrower V and 
smaller Ar/ d , as shown in Fig.8, where three models fi b = 0.02, 0.044, 0.09 are given. Therefore, 
a higher fi b leads to a lower Cf c , as shown in Fig.9 for these three values of fi b . 

The location of the peaks: 

From the formula (76), we also can analyze the peak's location of the power spectrums. 
The factor functions Pai[k(rjd — Vo)] in (77) and P^k^d — rjo)] in (78) are all combination of 
the spherical Bessel function ji(k(i]d — i]q)), which is peaked at I ~ k(i]Q — 7] d ) — kr]o for I 1. 
In Fig. 10 and Fig. 11 Pqi and Pci are plotted with I = 100, where it is shown that Pqi peaks 
at ki] ~ /, and Pa peaks at ki] ~ 1.27/. So the peak location of the power spectrums are 
directly determined by 
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C xx oc h( Vd ) k 2 D\k) \ k=l/vo , (82) 

The factor D(k) has a larger damping at larger /, so the first peak of the power spectrum 
has the highest amplitude. Let us look at the first peak of C xx , where D(k) ~ 1. Eq.(37) 
gives h{ri d ) 2 = Alk 2 (r] e /r] d ) 2 [A 1 j2(kr] d )+A 2 y2(kr] d )] 2 . Since j 2 (krj d ) term is the increasing mode 
and the y2{ki] d ) term is the damping mode for the waves inside the horizon, so we can take 
h(Vd) 2 ~ Alk 2 (r] e /rjd) 2 [Aij2(kr]d)] 2 . As j2(x) peaks at nearly x ~ 3, so h(r] d ) 2 peaks at kr] d ~ 3. 
Thus C xx peaks around 

l~kr] ~ 3r] /r] d . (83) 

This sudden-change approximation is a reasonably good estimate. If we use the WKB approx- 
imation in the ACDM Universe with r] d /r]Q = 0.0195, then we find that the factor function 
h{Vd) 2 k 2 m Eq.(82) is peaked around ki] Q ~ 127, shown in Fig. 12, while the sudden-change 
approximation in Eq.(83) gives kf]o ~ 154. Therefore, the estimate (83) holds approximately. 
The value of rjo/rjd is related to the dark energy component and baryon component. For in- 
stance, we take the three models with fi A = 0.65, 0.73, and 0.80, respectively, and with fixed 
fit = 0.044, Qdm = 1 — fiA — fib- Then a numerical calculation yields that r] /r]d — 50.1, 51.3, 
and 53.6, respectively. The gravitational waves h(r]d) also depends on fiA, as shown in Fig. 13, 
where it is clearly shown that a smaller Q\ will shift the peaks h(r]d) slightly to larger scales. 
Correspondingly, a smaller fi A will shift the peak of Cf x to larger scales, as demonstrated in 
Fig. 14. So we have the conclusion that a lower dark energy fi A makes the peak of C xx to 
locate at smaller /. This suggests a new way to study the cosmic dark energy. 

The baryon component also influence the decoupling time i] d . A larger fife has a larger 
decoupling time r] d , and therefore, a smaller / ~ 3r] /ri d . For fixed Q\ = 0.73 and fldm = 
1 — fiA — fife, the three models with fi b = 0.02, 0.044, and 0.09, respectively, are given in 
Fig.8. The corresponding values are r] /r]d = 54.9, 51.3, and 50.1, respectively. So we have 
the conclusion that a higher baryon density fife makes the peak to locate at smaller /, as is 
demonstrated in Fig.9. 

The influence of the spectrum index rir of relic GW on C xx : 

The initial condition in (32) of the the relic GW will influence the CMB polarization. Not 
only the value of initial amplitude h{0) but also the index rix in (33) will determine C xx . To 
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reveal how the index tit changes the Cf c , we plot in Fig. 15 three curves of Cf c for n T = —0.1, 
0.0, and 0.1, respectively, where the parameters are taken: r = 1, a = 2 in the half-gaussian 
model. Fig. 15 shows clearly that a larger value of ut produces a higher polarization spectrum 
(jcc fj^g reason f or this feature is the following: Cf c oc \h(i] d )\ 2 , and \h(i] d )\ differs from 
only at larger kr] > 50, as in Fig. 2, and \h(r} d )\ oc k\h(rj d )\. Since k 3 \h(r) d )\ 2 oc k nT , so in 
the range of larger k a larger value n-r will give a larger \h{rj d )\ 2 and a larger Cf c . Similar 
behavior also occurs for the electric polarization Cf G . 



8. Conclusion. 

In this paper we have analytically calculated the polarization power spectrum of CMB 
generated by relic gravitational waves in the Polnarev's method. As an approximate analytic 
result, it is quite close to the numerical one evaluated from the cmbfast code, especially, for 
the first three peaks of the spectrum that are observable. We have arrived at the analytic 
polarization spectra C xx with several improvements over the previous results. 

1. The spectrum C xx is proportional to the relic gravitational waves \h(r] d )\ 2 , the value of 
which is taken from the WKB approximation in our treatment, which is good compared with 
the numerical results. When looking for the location of the peaks of C xx analytically, we find 
that it is very convenient to employ the analytic exact solution \h(r] d )\ 2 from the sudden-change 
approximation. For instance, the first peak is found to be located at I ~ 3r] /r] d . Moreover, 
it has the merit that the dependence on the dark energy Q\ and the baryon Q b , through the 
decoupling time r] d , is also clearly demonstrated. Both lower dark energy and higher baryon 
density fib make the peak of C xx to locate at smaller I. 

2. The visibility function describing the decoupling process has been given by the analytic 
half-gaussian approximation, which fits analytically the actual decoupling process better than 
the usual Gaussian model. For example, the improvement of V{rf) on the left side of t] d in the 
half-gaussian model is considerable, which is about ~ 11.5%, and correspondingly, this also 
leads an improvement on the spectrum AC XX /C xx ~ 30% around the second and the third 
peaks. The time integration V(i]) yields the damping factor D(k), which contains a parameter 
a in the range [0, 2], coming from the phase of the product of the h{rj) and ji(k(i] — i] )). In 
particular, our half-gaussian model with a G (1.7,2.0) gives a reasonably good fitting to the 
first three peaks of the spectra C xx . 

3. In dealing with the Boltzmann equations analytically, we have worked up to the second 
order of the tight coupling limit, resulting in an amplitude of C xx smaller than that in the 
usual first order by 58%. 

5. We have studied the influence of the power spectrum index of the relic GW on the CMB 
polarization, and have found that a larger value of tit produces higher polarization spectra 
Cf c and Cf G . 
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Figure 1: The gravitational waves h(rjd) depends on the scale factor a{rf). The solid line is the result of the 
sudden-change approximation, the dash line is the WKB approximation result, and the dot line is the numerical 
results, which is nearly overlapped with the dash line. Note that, here the initial gravitational waves have been 
rescaled with h(k) = 1 for demonstrational purpose. 
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Figure 2: The gravitational waves h(i]d) depends on the scale factor 0(77). The solid line is the result of the 
sudden-change approximation, the dash line is the WKB approximation result, and the dot line is the numerical 
results, which is nearly overlapped with the dash line. Again the rescaled h(k) — 1 is used as in Fig.l. 
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Figure 3: The visibility function V(rj) around the decoupling. The solid line denotes the numerical result, the 
dash dot line denotes the gaussian fitting in Eq.(45) with Arjd/ilo — 0.00143 and V{r]d)r}o — 279, and the dot 
line denotes our half-gaussian fitting in Eqs.(46) and (47) with Ar]di/i]Q = 0.00110 and Ary^/jjo = 0.00176. 
Here r/d/r/o = 0.020. There is an area of strip on the left side of rjd in the Gaussian model differing from 
the numerical result, which is about ~ 11.5%. This variation will lead to an error in the spectrum Cf x 
correspondingly. 
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Figure 4: The "electric" polarization power spectrum C GG with the ratio r = 1. The solid line denotes 
the numerical spectrum from CMBFAST code, and the dot lines denote the analytic results. The upper dot 
line is the result of the half-gaussian visibility function with a = 1.7, the middle dot line is the result of the 
half-gaussian visibility function with a = 2 and, and the lower dot line is the result from the Gaussian visibility 
function with a = 2. While at large scales these models are close to each other, the half-gaussian model 
improves the spectrum height of the 3 rd peak by about A log 10 C GG ~ 0.2 . 
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Figure 5: The "magnetic" polarization power spectrum Cf c with the ratio r = 1. The solid line denotes 
the numerical spectrum from CMBFAST code, and the dot lines denote the analytic results. The upper dot 
line is the result from the half-gaussian visibility function with a = 1.7, the middle dot line is the result from 
the half-gaussian visibility function with a = 2, and the lower one is the result from the Gaussian visibility 
function with a = 2. Again at large scales these models are close to each other, but the half-gaussian model 
improves the spectrum height of the 2 nd peak by about Alog 10 Cf c ~ 0.1 . 
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Figure 6: The WMAP and Planck satellite measurements on the CMB magnetic polarization signal. The three 
solid curves show the analytic polarization power spectrums Cf c for the tensor-scalar ratio r — 0.3, 0.1, 0.01, 
respectively, in the ACDM Universe with Q b = 0.044, tt dm = 0.226, fi A = 0.73 . 




Figure 7: The WMAP and Planck satellite measurements on the CMB magnetic polarization signal. The 
parameters are the same as in Fig. 6. The solid lines show the numerical Cf c using cmbfast. Here we have 
included the influence of cosmic reionization with the reionization optical depth n r = 0.09. 
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Figure 8: The dependence of the visibility function V(rj) on the baryon Oh in ACDM Universe with Qa = 0.73, 
and Qdm = 1 — £1a — The baryon density has been taken to be fib = 0.02, 0.044, 0.09, respectively. 




Figure 9: The dependence of the magnetic polarization power spectrums on the baryon in ACDM Universe 
with Oa = 0.73, fldm — 1 — flA — Qb, an d r = 1. The baryon Q b = 0.02, 0.044, 0.09, has been taken respectively, 
as in Fig. 8. 
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Figure 10: The factor -P<|;(fc) as given in Eq.(77) with fixed I = 100 as a function of the wave number k. 
Obviously it is peaked around krjo ~ 100, verifying the relation (83). 




Figure 11: The factor P^j(fc) as given in Eq.(78) with fixed I — 100 as a function of the wave number k. It is 
approximately peaked around krj ~ 127, thus (83) as an estimate holds approximately. 
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Figure 12: In the WKB approximation the factor function \h(rjd)\ 2 k 2 rio/AT is peaked around <~ 127, validating 
the relation (83) as a fairly good estimate. Here the factor \h{rjd)\ 2 k 2 is multiplied by a factor tjo/At only for 
a clear graphical demonstration. 
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Figure 13: The influence of the dark energy Sl\ on the time derivative of GW, h{rjd), at the decoupling. A 
larger value of f^A shifts the peaks of h{rjd) to smaller scales. 
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Figure 14: The magnetic polarization spectrum Cf c depends weakly on Q\, and a larger value of Oa shifts 
the peaks of slightly to smaller scales. 
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Figure 15: The magnetic polarization spectrum Cf depends on the spectrum index tit of the relic GW. 
Cf' c are plotted for three values of nx = —0.1, 0.0, and 0.1, respectively. The following parameters are taken: 
r = 1, a = 2 in the half-gaussian model. A larger value of nr produces a higher spectrum Cf c . 
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